
% Idiosyncratic Volatility, Growth Options, and the Cross-Section of Asset
% Returns, by A. Barinov and G. Chabakauri
% 
% Model solution using Monte Carlo simulations.
% g and h are idiosyncratic and systematic volatility parameter vectors.
% Pictures are two-dimensional as functions of g and h,
% betas are also expressed as functions of g and P/D by running a regression of
% betas on the latter variables

clear;
clc;

% Exogenous Parameters
Gamma=1;
T=1; dt=0.000001; N=10; % N^2 is the number of stocks; 
K=50000; % number of simulations
r=0; %interest rate

e=ones(1,N);


% Calibration 1 (main calibration)
rho=0.8; 
v1_bar=0.4;
kappa1=0.05;
c1=-0.1;
v2_bar=0.4;
kappa2=0.05;
c2=-0.1;
mu_D=0.017;
Strike=1.2; 
lambda=2.5; % convexity
p=0.7; % weight on non-option tree, 1-p is weight on the option.

% Calibrations 2-5 are presented in robustness appendix

% Calibration 2 (in the robustness appendix)
%  rho=0.8;
%  v1_bar=0.3;
%  kappa1=0.02;
%  c1=-0.1;
%  v2_bar=0.3;
%  kappa2=0.02;
%  c2=-0.1;
%  mu_D=0.017;
%  Strike=1.2; 
%  lambda=2; 
%  p=0.7; 

% Calibration 3 (in the robustness appendix)
% rho=0.8;
% v1_bar=0.3;
% kappa1=0.05;
% c1=-0.1;
% v2_bar=0.3;
% kappa2=0.05;
% c2=-0.1;
% mu_D=0.017;
% Strike=1.2; 
% lambda=2.5; 
% p=0.7; 

% Calibration 4 (in the robustness appendix)
%rho=0.8;
%v1_bar=0.4;
%kappa1=0.05;
%c1=-0.2;
%v2_bar=0.4;
%kappa2=0.05;
%c2=-0.2;
%mu_D=0.017;
%Strike=1.2; 
%lambda=2.5; 
%p=0.7; 

% Calibration 5
% rho=0; 
% v1_bar=0.3;
% kappa1=0.05;
% c1=-0.1;
% v2_bar=0.3;
% kappa2=0.05;
% c2=-0.1;
% mu_D=0.017;
% Strike=1.2; 
% lambda=2.5; 
% p=0.7; 


%--------------------------------------------------------------------------

% Generating systematic and idiosyncratic volatility parameters
h=ones(1,N);
h0=1;
g=zeros(1,N);
for i=1:1:N
    g(i)=(i)*1.5/N;
end
for i=1:1:N
    h(i)=h0+0.1*(i-1)*1/N;
end


H=reshape(h'*e, N^2,1)';
G=reshape(e'*g, N^2,1)';

% Computing alphas, betas, and other parameters:
[v1_mean, v2_mean, alpha, alpha_hat_fm, beta_M_fm, beta_v1_fm, beta_v2_fm, S, SM, ER, ERM, Data]=Betas_2D(h0, H, G, Gamma, T, dt, r, N, v1_bar, kappa1, c1, v2_bar, kappa2, c2, mu_D, Strike, lambda, p, K, rho);

%==========================================================================
% Cross-sectional regressions of CAPM and factor alphas and IVol and systematic volatility betas on g and P/D

beta_v1_vec=reshape(beta_v1_fm, N^2,1)';
beta_v2_vec=reshape(beta_v2_fm, N^2,1)';
alpha_CAPM_vec=reshape(alpha, N^2,1)';
alpha_fm_vec=reshape(alpha_hat_fm, N^2,1)';
S_vec=reshape(S, N^2,1)';

% Range for P/D ratios
p=zeros(1,N);
p_min=min(S_vec);
p_max=max(S_vec);
for i=1:1:N
    p(i)=p_min+(i)*(p_max-p_min)/N;
end
P_matrix=p'*e;
G_matrix=e'*g;

X1=1+0.*G';
X2=G';
X3=S_vec';

X=[X1 X2 X3];
Omega=(X'*X)\X';
%--------------------------------------------------------------------------
% Beta_sys as function of g and P/D

a=Omega*beta_v1_vec';
R_sq_sys=std(X*a)/std(beta_v1_vec');
beta_sys_fitted=reshape(X*a, N, N);
beta_sys_reg=a(1)+G_matrix.*a(2)+P_matrix.*a(3);

% Beta_IVol as function of g and P/D
b=Omega*beta_v2_vec';
R_sq_ivol=std(X*b)/std(beta_v2_vec');
beta_IVol_fitted=reshape(X*b, N, N);
beta_IVol_reg=b(1)+G_matrix.*b(2)+P_matrix.*b(3);

%--------------------------------------------------------------------------
% alpha_CAPM as function of g and P/D

c=Omega*alpha_CAPM_vec';
R_sq_CAPM=std(X*c)/std(alpha_CAPM_vec');
alpha_CAPM_fitted=reshape(X*c, N, N);
alpha_CAPM_reg=c(1)+G_matrix.*c(2)+P_matrix.*c(3);

d=Omega*alpha_fm_vec';
R_sq_fm=std(X*d)/std(alpha_fm_vec');
alpha_fm_fitted=reshape(X*d, N, N);
alpha_fm_reg=d(1)+G_matrix.*d(2)+P_matrix.*d(3);

%==========================================================================

% Reporting results

surf(g,h, alpha);
title('CAPM alpha','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$h$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,h, alpha_hat_fm);
title('Factor alpha','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$h$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,h, beta_M_fm)
title('Market beta','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$h$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,h, beta_v1_fm)
title('Systematic vol beta','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$h$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,h, beta_v2_fm)
title('Idiosyncratic vol beta','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$h$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,h, S)
title('S/D ratio','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$h$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,p, beta_sys_reg);
title('Systematic vol beta','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$S/D$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,p, beta_IVol_reg);
title('IVol beta','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$S/D$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,p, alpha_CAPM_reg);
title('CAPM alpha','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$S/D$'},'fontsize',24,'interpreter','latex');
figure;

surf(g,p, alpha_fm_reg);
title('Factor alpha','fontsize',20,'interpreter','latex');
xlabel({'$g$'},'fontsize',24,'interpreter','latex');
ylabel({'$S/D$'},'fontsize',24,'interpreter','latex');

